Thermal Instability with Anisotropic Thermal Conduction and Adiabatic 
Cosmic Rays: Implications for Cold Filaments in Galaxy Clusters 



Prateek Sharma^, Ian J. Parrish^, Eliot Quataert 

Astronomy Department and Theoretical Astrophysics Center, University of California, Berkeley, 

CA 94720 

psharmaSastro . berkeley . edu , iparrishOastro . berkeley . edu , eliotSastro . berkeley . edu 



ABSTRACT 

Observations of the cores of nearby galaxy clusters show Ha and molecular emission 
line filaments. We argue that these are the result of local thermal instability in a 
globally stable galaxy cluster core. We present local, high resolution, two-dimensional 
magnetohydrodynamic simulations of thermal instability for conditions appropriate to 
the intracluster medium (ICM); the simulations include anisotropic thermal conduction 
along magnetic field lines and adiabatic cosmic rays. Thermal conduction suppresses 
thermal instability along magnetic field lines on scales smaller than the Field length 
(>10 kpc for the hot, diffuse ICM). We show that the Field length in the cold medium 
must be resolved both along and perpendicular to the magnetic field in order to obtain 
numerically converged results. Because of negligible conduction perpendicular to the 
magnetic field, thermal instability leads to fine scale structure in the perpendicular 
direction. Filaments of cold gas along magnetic field lines are thus a natural consequence 
of thermal instability with anisotropic thermal conduction. This is true even in the 
fully nonlinear regime and even for dynamically weak magnetic fields. The filamentary 
structure in the cold gas is also imprinted on the diffuse X-ray emitting plasma in 
the neighboring hot ICM. Nonlinearly, filaments of cold (~ 10^ K) gas should have 
lengths (along the magnetic field) comparable to the Field length in the cold medium 
~ 10~^ pc! Observations show, however, that the atomic filaments in clusters are far 
more extended, ~ 10 kpc. Cosmic ray pressure support (or a small scale turbulent 
magnetic pressure) may resolve this discrepancy: even a small cosmic ray pressure in 
the diffuse ICM, ~ 10~^ of the thermal pressure, can be adiabatically compressed to 
provide significant pressure support in cold filaments. This is qualitatively consistent 
with the large population of cosmic rays invoked to explain the atomic and molecular 
line ratios observed in filaments. 
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Introduction 



The thermal instabihty has been studied extensively in the context of the interstellar me dium 



(ISM; Field Ill965 : Koyama Sz Inutsuka 2000 : Sanchez-Salcedo. Vazquez-Semadeni. Sz Gazol ]Eo02; 



2OO5I ) and the formation 



Kritsuk &: Norman 



2002; 



Piontek Sz Qstriker 



2004 : 



Audit fc Hennebelle 



of solar prominences (e.g.. iKarpen et al. Ill988l ). but its role in galaxy clusters has not received as 
much attention. The cooling time in the intracluster medium (ICM) near the centers of galaxy 
clusters may be as s hort as 10-100 Myr. O bservations show that there is a dramatic lack of plasma 
below ~ 1 keV (e .g.. iPeterson et al. 1120031 ). inconsistent with the prediction of the original cooling 



flow models (e.g.. lFabian I ll994l 'l. In addition, the star formation rate in the ceiitral g alaxy is 10-100 



times smaller than if the gas cooled at the predicted rate (e.g., lO'Dea et al. II2008I ). This implies 
that cooling in the intracluster medium (ICM) is balanced by some form of heating that maintains 
an approximate global thermal equilibrium. Feedback from a central Active Galactic Nucleus 
(AGN ) is an energetically plausible source of the required heating (e.g., iGuo. Oh. &: Ruszkowski 



20081 ). However, precisely how the AGN provides this heating is not understood in detail; nor are 
other heating mechanisms ruled out. Many clusters with a short central cooling time (< 1 Gyr; or 
equivalently a low centr al entropy) show both star formation and Ha emission, indicative of cool 
plasma at < 10^ K (e.g.. lCavagnolo et al. II2OO8I ). Even if heating balances cooling in a global sense, 
the ICM plas ma is expec ted to be locally thermally unstable because of the form of the cooling 
function (e.g.. iField Ill965l ). This local thermal instability is an attractive mechanism for producing 
the Hq and molecular filaments seen in clusters with short cooling times. 

Recently, the ato mic and molecular filaments in the core of the Perseus cluster have been 



spatially resolved (e.g.. IConselice et al. 



coherence of these filaments. 



2001 



Fabian et al. 



horn 



Salome et al. 1120061 ). Based on the narrowness and 



suggested that magnetic fields play a critical role 



in the dynamics of the filaments. In addition to the possible role of magnetic pressure and tension, 
the magnetic field also modifies the microscopic transport processes in the ICM because the mean 
free path along magnetic field lines is orders of magnitude larger than the gy roradius. As a result, 
thermal conduction is primarily along magnetic field lines (jBraginskii Ill965l ). 



This paper centers on carrying out magnetohydrodynamic (MHD) simulations of thermal in- 
stability with thermal conduction along magnetic field lines. We focus on understanding the physics 
of the thermal instability in the ICM, rather than on making detailed comparisons with observa- 
tions. Throughout this paper we ignore the possible presence of a background gravitational field. 
This allows us to study the physics of the thermal instability without the added complication of 
buoyant motions, infiow, etc. In addition to including the effects of anisotropic heat transport, 
we also include cosmic rays as a second fluid. Even an initially small cosmic ray pressure can be- 
come energetically important in cold filaments. Part of the motivation for including cosmic rays is 
that a significant population of energetic ions (» eV, the temperature of o ptical filaments) app ear 
required to explain the observed line ratios in filaments in galaxy clusters (jFerland et al. II2009I ). 



The multiphase nature of the ICM is physically analogous to the well-studied multiphase ISM. 
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The ISM has three dominant phases: a molecular phase at ~ 100 K, an atomic phase at ~ 10^ K, and 
the hot phase at 10® K. The cooling function in the ISM is thermally bistable with thermally stable 
phases at ~ 100 K and ~ 10^ K. The hot phase is thermally uns table but is probably maintained at 



its temperature by supernova heating (jMcKee &: Ostriker 1119771 ) . The same physical considerations 



apply for the ICM, except that the hot phase is maintained by a still poorly understood heating 
process (e.g., AGN feedback). 

This paper is organized as follows. Section 2 summarizes our model equations and the results of 
a linear stability analysis including conduction along field lines, cosmic rays, and magnetic fields (see 
appendix A). Section 3 presents the numerical set-up and the results of our numerical simulations. 
Section 4 discusses the astrophysical implications of our results. 



2. Governing Equations and Numerical Methods 

A magnetized plasma with cosmic rays can be described by the following two-fluid equations: 

dv , ^2 (B-V)B 

p_ = _V(, + ,„+_) + L_^, ,2) 

— = Vx{vxB), (3) 

§ - - ^ = -nen^A{T) -V-Q + H{t), (4) 
dt p dt 

and 

dt p dt ' ' ^ ^ 

where d/ dt = d / dt + v ■ V is the Lagrangian time derivative, A(T) is the cooling function, 

Q = -K||b(b- V)r (6) 

is the heat flux along magnetic fleld lines, 

T = -D\\b{b-V)p,, (7) 

is the diffusive cosmic-ray energy flux (multiplied by [7cr — 1] ) j p is the mass density, rig and rii are the 
electron and ion number densities respectively, v is the common bulk-flow velocity of the thermal 
plasma and cosmic rays, B is the magnetic field, b = B/B, p and pcr are the thermal-plasma 
and cosmic-ray pressures, is the parallel thermal conductivity, is the diffusion coefficient 
for cosmic-ray transport along the magnetic field, and 7 = 5/3 and jcr = 4/3 are the adiabatic 
indices of the thermal plasma and cosmic rays, respectively. We assume one-third solar metallicity 
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so that the mean molecular weights are fx = 0.62 and /ig = 1-18. We do not include gravity in the 
momentum equation (Equation ^) in order to focus on the thermal physics. Note that our model 
equations also do not include the streaming of cosmic rays re lative to the thermal plasma, which 



provides a mechanism for heating the thermal plasma (e.g., iLoewenstein. Zweibel. &: Begelman 



1991I : IGuo fc Oh 1 120081 ). This is numerically subtle to include (jSharma. Colella. k Martin I l2009l ) 



and will be studied in future work. 

It is difficult to study the problem of thermal instability without a well-defined equilibrium 
state. In order to ensure that we have such a state, at each time step the heating term H{t) in 
equation (jj]) is updated so that the volume averaged heating and cooling in our computational 
domain balance each other. Without such heating, the plasma as a whole cools to very low tem- 
peratures on a cooling time (the same timescale on which the thermal instability is developing). 
Since the source of heating and its functional form are not that well- understood in the ICM, we 
choose a constant heating per unit volume for simplicity. Calculations with a constant heating per 
unit mass, i.e., H{t) oc p, yield very similar results because cooling (oc n^) dominates in the cold 
phase and heating dominates in the hot phase in both cases. 



2.1. Linear Stability 

In appendix A we study the linear thermal stability of a uniform plasma with magnetic fields, 
cosmic rays, and thermal conduction along magnetic field lines. To isolate the physics of interest 
in this paper, we focus on the "condensation mode," i.e., t he entropy mode. This calculation 



is a straightforward generalization of previous results (e.g., iField I Il965l ). but we include it for 
completeness. Here we quote the final results. 

When the cosmic ray and magnetic pressure are negligible compared to the plasma pressure, 
and the cooling time tcooi is long compared to the sound-crossing time, the growth rate for the 
thermal instability is given by 

7 = -X||/c|-C^^, (8) 

where t~^^^ = (7 — l)neniA/p. Note that the thermal conductivity, k\\, in Equation ^ is related to 
the diffusivity used here, X||; by Ky = rieksXn- The first term on the right hand side of Equation 
([8]) describes the conductive stabilization of modes with short wavelengths parallel to the local 
magnetic field (large A;||). This implies that the fastest growing modes will be elongated along the 
magnetic field lines and hence filamentary. The critical parallel length-scale at which 7 = (the 
Field length) is given by: 



Xll^cool 



(iln(r2/A)/dlnT 



(9) 



When the cosmic ray and/or magnetic pressure is large compared to the plasma pressure, the 
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isochoric growth rate applies, i.e., 



9 1 din A , ^ 

7=-X„^f-Coi^- (10) 



This is also applicable when the cooling time is shorter than the sound-crossing time, irrespective 
of the magnetic and cosmic-ray contributions to the total pressure. However, for typical ICM 
conditions, the cooling time in the hot plasma is longer than the sound-crossing time. 



2.2. Thermal Conductivity and the Cooling Function 



The thermal conductivity of a fully ionized plasma is governed by electron collisions with the 
background ions and electrons. We are interested in plasmas hotter than 10^ K, so that 



1.84 X 10" 
In A 



_ 2^5/2 



erg s 



lK-7/2 



cm 



(11) 



where we use a Coulomb logarithm of In A = 37. Since the Larmor radius of electrons is much 
smaller than their collisional mean free path, thermal conduction is only effective along magnetic 
field lines; the perpendicular transport is negligible. 



We use the cooling function for an ionized plasma from [Sutherland &i Dopita I (jl993l ). A fit to 
th e Sutherland and Dopit a cooling rates for a third solar metallicity is given by (a generalization 
of ItozzI fc Norman |[200ll '): 



A 



10 



^2 (8.6 X lO^^Tj^^y + OmST^^^ + 0.063) ergss"^ cm^ 



keV 

for T > 0.02 keV, 
A = 6.72 X 10^22 (Tkcv/0.02)°-'' ergs s-^cm^ 

for T < 0.02 keV, T > 0.0017235 keV, 
A = 1.544 X 10'22 (Tkcv/0.0017235)^ ergs s^^cm^ 

for T < 0.0017235 keV, 



(12) 



where T^eV is the temperature in keV. For the cooling function given by Equation (I12p . the only 
phase that is thermally stable according to Equation ([8]) is plasma with T < 0.0017 keV ~ 10^ K. 



Kovama &: Inutsuka I (|2004l ) have shown that with isotropic conduction, thermal instability 



simulations do not converge with increasing resolution unless the Field length is always resolved by 
a few grid cells. We find the same result in our simulations. The Field length can be written as 



Ai. « 14.4 T^/^ (ne,o.ini,o.i)"'/' A„^/' 



dln(rVA) 
dlnT 



-1/2 



kpc 



(13) 



where ne^o.i (^i,o.i) is the electron (ion) number density in units of 0.1 cm and A = 10 A_23 erg 
cm^. Because the Field length decreases rapidly with decreasing temperature, it is prohibitive 
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to resolve the Field length in a numerical simulation if the plasma has a wide range of temperatures. 
Indeed, at fixed pressure, Equation (jl3p implies that the Field length at a few 10^ K is ~ 10*^ times 
smaller than at 1 keV! In order to ensure that our simulations always resolve the properties of the 
cold phase of the ICM, we use an artificial cooling function in which the thermally stable phase is at 
a much higher temperature of 2 x 10^ K. Because our simulations begin with plasma at typical ICM 
temperatures 10'' K, the modest range of temperatures on the computational domain makes it 
feasible to always resolve the Field length. Figured] shows a comparison of the true cooling function 
from Equation (|12|) (solid line) and our modified cooling function (dotted line). 

With anisotropic thermal conduction, we find that numerical convergence requires resolving 
both the parallel and perpendicular Field lengths in the cold stable phase. More precisely, if the 
perpendicular Field length is not resolved, the width of cold structures perpendicular to the local 
magnetic field decreases with increasing resolution. In order to ensure that there are no spurious 
results due to unresolved structures at the grid scale, we thus include a constant isotropic diff^usivity 
of 3 X 10^^ cm^s~^, in addition to the parallel thermal conductivity given by Equation (jlip . The 
perpendicular diff^usivity required for numerical convergence depends on resolution; we have verified 
numerically that ~ 3x 10^^ cm^s~^ is the minimum isotropic diffusivity required to obtain converged 
results for our two-dimensional simulations presented in ^ The thermal diffusivity along magnetic 
field lines (Equation (|lip ) is ~ 30 times larger than the isotropic diffusivity for our typical initial 
conditions. Thus thermal conduction is still primarily along the magnetic field, although the 
perpendicular conductivity is orders of magnitude larger than the microscopic value. In ^4.31 we 
discuss how this artificially large perpendicular conductivity might affect the conclusions drawn 
from our simulations. 



2.3. Numerical Simulations 



In this paper, we carry out local one and two dimensional numerical simulations of thermal 
instability for conditions appropriate to the ICM. We use unstratified local patches of the ICM 
to isolate the physics of the therm al instability, as opposed to th e buoyancy instabilities present 
in a stratified, conductive plasma (jBalbus Il2000l : iQuataert II2008I ). We focus on two-dimensional 
simulations - as opposed to three-dimensional simulations - because of the challenging numerical 
requirement of resolving the Field length in the cold medium (both parallel and perpendicular to 
magnetic field lines). Our box size is 40 kpc, somewhat larger than the typical Field length in the 
hot medium; we use periodic boundary conditions. 



We use the pubhcly available ZEUS-MP code (jHaves et al. 1 120061 ) to solve the MHD equa- 
tions. Thermal conduction along magnetic field lines is treated explicitly, using the method of 
Sharma &: Hammett I (|2007l ) to prevent unphysical negative temperatures. Since the stable timestep 
for conduction is much smaller than the MHD timestep, thermal conduction is subcycled. The 
cooling and heating terms in Equation ([4]) are combined at each grid point, and internal energy 
is updated by a first order explicit (semi-implicit) method if heating (cooling) dominates; this en- 



-7- 



sures that the internal energy is always positive, irrespective of the timestep. Since the typical 
cooling time is much longer than the sound-crossing time across a grid cell, this first order accurate 
treatment of cooling is sufficient. 

Our initial condition consists of plasma with T = 0.78 keV and ng = 0.1 cm~^; these parameters 
are characteristic of a reasonably dense, low entropy (~ 3.6 keV cm^) part of the ICM at small 
radii, deep in the cluster core. Note that the cooling time is longer than the sound crossing 
time across the box so that the thermal instability is in the roughly isobaric limit (except in 
magnetic/cosmic ray dominated regions, where it behaves isochorically; compare Equations (jlOp 
& dH])). We initialize homogeneous and isotropic ^ 1% isobaric density /temperature perturbations 
on this initial equilibrium state; the spectrum of initial perturbations is oc A; for k < kQ and oc k^^ 
for k > ko, so that most of the power is initially at ~ ko; /cq corresponds to a scale In/ko 0.8 
kpc for most of the simulations (we also experimented with smaller and larger ko for comparison). 
The exact spectrum of initial perturbations is somewhat arbitrary and is not well constrained in 
the ICM, but is also not crucial for the subsequent evolution. All of the simulations with the 
same kg have identical initial conditions and the initial conditions are smoothed/interpolated for 
lower /higher resolution simulations. This is important for quantitative testing of convergence with 
increasing resolution. The one-dimensional simulations do not include magnetic fields, while the 
initial magnetic field is B = 5 /xG and aligned at 45*^ to the box in the two-dimensional simulations. 
The simulations with cosmic rays begin with a uniform cosmic ray pressure. 

Tables [1] and [2] summarize the one and two dimensional simulations discussed in detail in this 
paper. The tables list the mass {fm) and volume (/y) fractions of plasma having a temperature 
below 5 X 10^ K in the nonlinear state; this is a reasonable proxy for the mass and volume fractions 
of the cold phase. To quantify the spatial coherence of the structures in the nonlinear state of the 
thermal instability, we define the parallel and perpendicular length scales of the density field via 

" f\b-V6p\dV 

and 

. M , (15) 

J\{z xb) -VdpldV 

where 6p = p — (p), {p) is the volume averaged density (which is constant in time because the mass 
in the computational domain is conserved), and z is perpendicular to the simulation plane. For 
the one-dimensional simulations x is used instead of b in Equation (I14p and L±_ is not defined, so 
we only provide -^^ii in Table [TJ 



3. One-dimensional Simulations 



Figure [2] shows temperature profiles in the linear {left) and nonlinear {right) regimes for one- 
dimensional hydrodynamic simulations with (HWC) and without (HNC) thermal conduction. For 
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Fig. 1. — Cooling function vs. temperature. Solid line: fit to the Sutherland & Dopita cooling 
function for a third solar metallicity (Equation (jl2p ). Dotted line: the modified cooling function 
used in our simulations which has a stable cold phase at < 2 x 10^ K. 
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Fig. 2. — Temperature profiles for one-dimensional runs with (HWC) and without (HNC) conduc- 
tion at different times in the linear {left) and nonlinear {right) regimes. The initial temperature 
fluctuations have been multiplied by a factor of 100 for clarity. 
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the simulation without conduction (dotted hne) the temperature fluctuations grow at all scales 
in the linear regime. By contrast, for the run with conduction (dashed line) modes with scales 
smaller than the Field length are suppressed by thermal conduction and only the large scale modes 
grow. Nonlinearly, the cold phase is compressed into a smaller and smaller volume with time in the 
absence of conduction, until the cold phase is unresolved (notice that in Table [J is much smaller 
for the run without conduction compared to the runs with conduction); the cold peaks also merge, 
reducing the total number of dense peaks in time (compare profiles at 1.43 and 1.91 Gyr for HNC 
in the right panel of Figure ED . Eventually all of the cold peaks will merge and approach the grid 
scale because there is no heating of the cold phase to prevent this. As the cold phase accumulates 
more and more mass in time, the negligible mass in the hot phase becomes hotter and hotter to 
conserve energy (which is enforced in our simulations via the heating term H{t) in Equation Q). 

The nonlinear evolution with conduction is qualitatively different: there is only one cold region 
(this is because of the large Field length in the initial plasma; this result is insensitive to the initial 
density fiuctuation spectrum) and the hot phase saturates at a temperature ^ 2 keV, much cooler 
than in the simulations without conduction. The temperature of the cold phase is, however, the 
same with and without conduction; this is set by the temperature of the thermally stable branch 
of the cooling function. Figure [2] shows that in the presence of conduction, the temperature profile 
reaches an approximate steady state, with very little change from 1.43 to 1.91 Gyr. The steady state 
requires both the additional heating H{t) in Equation ^ and thermal conduction. In particular, 
the cooling is dominated by the dense, cold gas while the extra heating is primarily supplied to the 
hot phase (because H is constant per unit volume). This extra heating is conducted to the rapidly 
cooling (oc n^) cold phase producing a steady state. This energy transfer from the hot to the cold 
phase can only be properly captured if the Field leng th in the cold phase is reso lved, which is why 



it is critical to do so to obtain converged results (see lKoyama &: Inutsuka II2004I ) 



Figure [3] shows temperature profiles at 1.43 Gyr for simulations including thermal conduction 
at several different resolutions. The temperature profile is reasonably converged only for simulations 
with more than 512 grid points; in particular, note that the physical size of the cold phase does 
not change with resolution for N > 512. This is only true when the Field length in the cold phase 
is resolved. The Field length for the initial temperature and density is ~ 10 kpc. The Field length 
for the isobaric cold phase at 2 x 10^ K (the stable phase for our modified cooling curve; Figure [1]) 
is ~ 0.07 kpc; this just starts to be resolved at more than 512 grid points since our box size is 40 
kpc, and thus Ax = 0.078 kpc at iV = 512. 



4. Two-dimensional Simulations 

Having used one-dimensional simulations to describe the basic physics of the thermal instability 
and the numerical requirements for simulating it, we now turn to the more physically realistic case of 
two-dimensional simulations. As we have emphasized previously, the Field length must be resolved 
both parallel and perpendicular to the local magnetic field in multi-dimensional simulations of the 



- 10 - 



Table 1: One dimensional runs 



LabeP 


Res. 


Ax (kpc) 


L||(kpc)'' 


J m 


fv 


HWC 


1024 


0.039 


3.26 


0.6 


0.043 


HNC 


1024 


0.039 


0.22 


0.9 


0.07 


HWCl 


512 


0.078 


3.26 


0.61 


0.043 


HWCh 


2048 


0.02 


3.26 


0.6 


0.042 


HWCll 


256 


0.16 


0.89 


0.67 


0.043 



"H stands for hydro. WC for with conduction. I & h stand for lower and higher resolution runs. Initially =0.1 
cm~^, T = 0.78 keV, so that the cooling time is ~ 95 Myr. The Field length in the initial condition is ~ 10 kpc and 
at 2 X 10*' K (temperature of the stable phase) is « 0.07 kpc. The box size is 40 kpc. 

''L|l is defined in Equation (|14p and is evaluated at 1.43 Gyr, when the results have reached a quasi-steady state. 
fm (fv) is the mass (volume) fraction of plasma below 5 x lO^K (the "cold phase") evaluated at 1.43 Gyr. 




Fig. 3. — Temperature profiles for one-dimensional simulations with conduction at t = 1.43 Gyr for 
different resolutions: HWCll (256); HWCl (512); HWC (1024); and HWCh (2048). Convergence is 
achieved for > 512 grid points. 
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thermal instability. This is why we (1) include an isotropic thermal conductivity (which helps 
resolve structures perpendicular to the field; see ^2.2p in addition to the parallel conductivity and 
(2) artificially increase the temperature of the thermally stable phase (Figure [1]). 

4.1. The Fiducial Run: MHD with Anisotropic Thermal Conduction 

Our fiducial two-dimensional simulation is MWC summarized in Table [21 this is an MHD 
simulation with anisotropic thermal conduction and an initial magnetic field of B = 5/iG aligned 
at 45" relative to the x-axis, in the plane of the simulation. The initial density (0.1 cm~'^) and 
temperature (0.78 keV) correspond to an initial cooling time of 95 Myr and a Field length ~ 10 
kpc along the magnetic field and ~ 1.7 kpc perpendicular to the field; the magnetic field initially 
contributes only ~ 0.4% (/? = Snp/B^ ~ 250) of the total pressure. Figure U] shows contour plots 
of the temperature in the linear (0.475 Gyr) and nonlinear (0.95, 1.425 Gyr) regimes, along with 
arrows showing the magnetic field direction at each time; because the pressure remains relatively 
constant even nonlinearly, density scales nearly as the inverse of temperature. 

Figure H] shows that the thermal instability develops anisotropically, with a filamentary struc- 
ture along the magnetic field; this is because thermal conduction efficiently suppresses small-scale 
structures along the field, but not across it. Quantitatively, the ratio Ly/L^ measures the anisotropy 
of filaments with respect to the magnetic field; this is ~ 2.5 at 0.95 Gyr (Table [2|) and increases 
to ~ 3.5 at later times (see the left panel in Figure [TTJ. The number of cold filaments decreases 
in time because some of the filaments merge together nonlinearly. Interestingly, the majority of 
the cold filaments are oriented along the direction of the local magnetic field even in the nonlinear 
regime. Some of the filaments at 0.95 Gyr are quite small and relatively isotropic because of small 
and nearly isotropic conduction in the cold phase. However, at later times (e.g., 1.425 Gyr) the 
small filaments coalesce to form large ones. The nonlinear development of the thermal instability 
proceeds in two phases: in the first phase nonlinear filaments aligned along field lines condense 
from the hot ICM, becoming shorter in time because of a smaller conductivity in the cold phase; 
in the second stage these cold filaments with large velocities (primarily along themselves) merge to 
form longer filaments. This is clearly seen in Figure [TT] as an increase in L^^/L± after an initial dip 
at ~ 1 Gyr. 

Figure H] shows that the direction of the magnetic field is only moderately perturbed from its 
initial direction even in the fully nonlinear regime. However, the magnetic field strength increases 
by a factor of > 3 — 8 in the cold filaments (see the left panel of Figure [5]), to the point where 
the magnetic pressure is important in the filaments. The regions over which the field is enhanced 
are coincident with, but significantly longer than, the location of the cold filaments. The field 
enhancement occurs via flux freezing as the cooling plasma is compressed perpendicular to the 
initial field direction in the nonlinear state of the thermal instability; analogous compression along 
the field lines is suppressed because of thermal conduction. In the hot diffuse gas between the 
filaments, the magnetic field decreases by a factor ~ 2 — 3 from its initial value of ~ 5 /xG. Note 
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Table 2: Two dimensional runs 



LabeP 


Res. 


Ax = Ay (kpc) 


Pct/P 


(cm^s-^) 


n T"- ^ — 

(heating) 

(cooling) 




J m 


fv 


MWC 


1024 


0.039 







1 


2.48 


0.51 


0.061 


MWCl 


512 


0.078 







1 


2.18 


0.5 


0.06 


MWCh 


2048 


0.02 







1 


2.55 


0.51 


0.064 


MWIC 


1024 


0.039 







1 


1.08 


0.48 


0.041 


MWCCR 


1024 


0.039 


0.1 





1 


3.45 


0.41 


0.12 


MWCCRs 


1024 


0.039 


10-3 





1 


2.5 


0.5 


0.063 


MWCCRd28 


1024 


0.039 


0.1 


1028 


1 


2.66 


0.4 


0.1 


MWCCRd30 


1024 


0.039 


0.1 


1030 


1 


2.29 


0.47 


0.061 


MWChO.gc 


1024 


0.039 







0.9 


3.91 


0.9 


0.38 


MWChl.05c 


1024 


0.039 







1.05 


3.86 


0.018 


1.5 X 10-3 



"M stands for MHD. WC means with conduction, IC is for isotropic conduction, CR for cosmic rays. 1 & h stand 
for lower and higher resolution runs. All runs have a small isotropic conduction added for convergence (see tj2.2|l . 
Initially Ue = 0.1 cm~^ and T = 0.78 keV, so that the cooling time is ~ 95 Myr. Initial magnetic field is 5 and 
aligned 45° to the two-dimensional cartesian box. The Field length in the initial condition is « 10 kpc and at 2 x 10® 
K (temperature of the stable phase) is « 0.07 kpc. The box size is 40 kpc. Some less crucial simulations are not 
included in the table but are discussed in the text. 

*The fiducial run. 

and L± are defined in Equations (|14p & (|15p . and are evaluated at 0.95 Gyr. 
'^fm ifv) is the mass (volume) fraction of plasma below 5 x lO^'K (the "cold phase") evaluated at 0.95 Gyr (except 
for MWIC where these are evaluated at 1.43 Gyr). 




x(kpc) x{kpc) x(kpc) 



Fig. 4. — Contour plots of Loqiq temperature (in keV) for the fiducial run (MWC) at linear (0.475 
Gyr; lef£) and nonlinear (0.95 Gyr, center, 1.425 Gyr, right) stages of the instability. The arrows 
show the magnetic field direction. 
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that for a realistic cooling function, the density contrast between the filaments and the diffuse 
medium will be larger than is found in our simulations, and so the magnetic field compression in 
the filaments will also be stronger. 

The right panel of Figure [5] shows that the velocities driven by the thermal instability can 
reach 30 — 100 km s~^, comparable to the sound speed in the cold filaments, but much less than 
the sound speed in the hot phase. Such high velocities can disrupt the tendency of buoyancy 



instabilities in the hot phase of the ICM to reorient the magnetic field (e.g.. ISharma et al. 112009b 



Parrish. Quataert. &: Sharma Il2010l ) . The high velocities are spatially coincident with the magnetic 
field enhancements and the cold filaments. The velocity vectors generally point toward the cold 
filaments in the hot phase, showing that mass from the hot thermally unstable medium is condensing 
into the cold phase. This fiow of mass is, however, transient. The thermal instability reaches a 
steady state in which cooling from the dense, cool ICM is balanced by conductive heating from 
the hot ICM, which is in turn heated (artificially) by our external heat source H{t) in Equation 
(jl]). Once this steady state is established, mass fiow between the phases is significantly reduced. 
Although mass flow across the phases is reduced, the cold filaments retain large velocities along 
themselves and the volume averaged velocity is ~ 20 kms~^ (see the right panel of Figure [11] 
discussed later). 

Nonlinearly, the plasma exists in two phases, with very little plasma at the intermediate 
temperatures. Figure [6] shows the mass {left panel) and volume (right panel) distribution of plasma 
at different times for the fiducial run. The plasma is at ~ 10^ K initially but evolves into a two- 
phase structure. The phase structure evolves rapidly at early times (before ~ 1 Gyr), but the 
evolution is slower at later times. The mass and volume occupied by the plasma at intermediate 
temperatures decreases in time. The "mass dropout rate," (i.e., the rate at which plasma cools 
below a given temperature) at 10^ K is large initially, but once a two-phase medium is established, 
the mass and volume of the hot and cold phases are roughly constant in time, with very little mass 
dropout. While there is significant mass in the cold filaments, most of the volume is occupied by 
the hot phase (see fm and fv in Table [2]) . The hottest plasma in the domain slowly becomes hotter 
with time in the two dimensional simulations; by contrast, in ID the plasma reaches a steady state 
at 1.43 Gyr (Figure [2]). It takes longer to reach a quasi-steady state in two dimensions because it 
is easier for hot isothermal regions to become thermally isolated from the cold plasma (because of 
the small perpendicular conductivity). Since the hottest plasma becomes hotter with time and the 
conductivity is a strong function of temperature (Equation ([TT]) ). it becomes difficult to run the 
simulations for long times. 



4.2. Simulations with Isotropic Thermal Conduction 

To assess the importance of including anisotropic thermal conduction, we carried out simu- 
lations identical to the fiducial run in every way except that the conductivity is isotropic at the 
Spitzer value (MWIC in Table [2]). Figure [7] shows the temperature contour plots at 0.475 Gyr {left 
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Log^^lBI(fxG) 



Log^Jvl(kms ^) 




20 
x(kpc) 



20 
x(kpc) 



Fig. 5. — Contour plots showing Logio\B\ (magnitude of the magnetic field strength) (left) and 
LogiQ\v\ (magnitude of the velocity) (right) for the fiducial run at 0.95 Gyr. The arrows in the 
velocity plot show the direction of the velocity unit vector. 
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Fig. 6. — The mass {dM / dlogi^T; left) and volume ( dV / dlogi^T; right) fractions occupied by 
plasma of a given temperature T for the fiducial run (MWC) at different times. The normalization 
is such that the total mass/volume under the curve is unity. The initial cooling time is ~ 0.1 Gyr 
and the simulations begin to saturate after ~ 0.8 Gyr. The hottest plasma in the box becomes 
hotter with time. 
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panel) and 0.95 Gyr {right panel). In the linear state the modes are isotropic and on relatively large 
scales, irrespective of the magnetic field direction. By contrast, with anisotropic conduction, the 
cold plasma is filamentary even in the linear state (Figure Nonlinearly, the orientation of the 
cold plasma in simulations with isotropic conduction is unrelated to - or even somewhat perpendic- 
ular to (see dotted line in Figure [8|) - the local magnetic field direction, unlike in simulations with 
anisotropic conduction, where the filaments develop along the magnetic field (Figured]). Although 
the morphology of the cold gas is different in the two cases, the evolution of the phase structure is 
qualitatively similar; there is significant mass in the cold phase, but the volume is dominated by the 
hot phase. The differences between Figures S] and [7] emphasize the critical importance of including 
anisotropic thermal conduction when studying the thermal physics of galaxy cluster plasmas. 



4.3. Convergence of Two-dimensional Simulations 

As described previously, in multi-dimensional simulations, the Field length must be resolved 
both along and perpendicular to the direction of the magnetic field in order for the numerical 
results to converge. Figure [9] shows temperature contour plots at 0.95 Gyr for runs including 
perpendicular conduction, with 2048 and 512 grid points, respectively. The temperature contour 
plots are reasonably similar, and are similar to the results for = 1024 in Figure HI Figure [TT] 
provides a more quantitative test of the convergence of the simulations: it shows that the volume 
averaged values of L||/L_|_ (the anisotropy of the filaments) and {\v\) (the random velocity) are 
almost the same, irrespective of resolution, for runs with perpendicular conduction (labeled "Y"). 

To explicitly illustrate the importance of including thermal conduction perpendicular to field 
lines for convergence, we carried out simulations similar to the fiducial run, but without the small 
isotropic conductivity. Figure fTOl shows temperature contour plots for simulations without perpen- 
dicular conduction, for A^ = 2048 and 512 grid points. In this case, the two different simulations 
give very different results; in particular the filaments are much thinner and the number of filaments 
is much larger (by a factor ~ 4) for the higher resolution run. This is also seen in the left panel of 
Figure [TT] which shows that the anisotropy of the filaments L||/L_|_ increases with increasing resolu- 
tion for runs without perpendicular conduction. Similarly, the volume averaged velocity ((|f |)) does 
not show clear convergence with increasing resolution in the absence of the isotropic conductivity 
(see the right panel of Figure [TT]) . 

It is important to stress that the simulations with perpendicular conduction (e.g.. Figure 
[9]) significantly over-estimate the thickness of the filaments perpendicular to the magnetic field. 



^The Field length perpendicular to the magnetic field is much smaller in the simulation with anisotropic conduction 
than in the simulation with isotropic conduction. This is why there is much more small-scale structure, and more 
cold 'filaments,' in Figure[3]than in Figure[71 In addition, because we initialize power primarily at « 0.8 kpc f tj2.3|) . 
the amplitude of the initial perturbations that can actually grow (> the Field length) is larger in the simulation with 
anisotropic conduction. These perturbations thus evolve somewhat more rapidly. 
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because the perpendicular conductivity is too large by orders of magnitude. In this sense the trend 
in Figures [11] and [10] is correct, namely the perpendicular structures should indeed be thinner than 
in our fiducial simulations. However, it is critical that the Field length be resolved perpendicular to 
the magnetic field, or else spurious numerical results can arise (e.g., in our simulations with cosmic 
rays, we found that the cosmic ray pressure could become spuriously large in cold filaments when 
they were not properly resolved, even if the initial cosmic ray pressure was negligible). Physically, 
small scale turbulent heat transport (e.g., due to Kelvin-Helmholz instabilities at the boundaries 
of the filaments) or other physics (e.g., cosmic-ray pressure; ^4.4p probably sets the perpendicular 
scale of the filaments, not the true microscopic perpendicular heat transport. These processes are 
not currently well-understood and it is unclear to what extent they can simply be treated as an 
enhanced perpendicular conductivity (as we have done here). 



4.4. Effects of Cosmic Rays on Filament Formation 

In the previous sections, we have highlighted the dynamics and thermodynamics of the thermal 
plasma during thermal instability. In this section we consider the role of cosmic-rays, i.e., a non- 
thermal population of particles. Figure [12] shows contour plots of the ratio of the cosmic ray to 
plasma pressure for simulations with two different initial cosmic ray pressures, Pcv/p = 0.1 and 10~^, 
respectively; the cosmic-rays are adiabatic in these simulations. Figure [12] shows that the cosmic 
rays become concentrated in the cold filaments; this is because the cosmic ray entropy Pcr/ P^^'^ is 
conserved and the cosmic rays are thus compressed along with the thermal plasma into the cold 
filaments. For the simulations with a very small cosmic-ray pressure [right panel of Figure [T2]) . the 
properties of the thermal plasma in the filaments and in the diffuse ICM are very similar to those 
in the simulations without cosmic rays (Figure [1]). In particular, because the cosmic ray pressure 
is small even in the nonlinear state, the cosmic rays do not affect the physics of how the filaments 
form. On the other hand, when the initial cosmic ray pressure is larger [left panel of Figure [T2]l . 
adiabatic compression of the cosmic rays in the filaments leads to cosmic ray pressure dominated 
filaments that are longer and broader than in the absence of cosmic rays; the additional cosmic ray 
pressure halts the contraction of the filaments when p^r ~ p- 

Table [2] shows that the volume fraction fy of the cold phase is larger for simulations in which the 
filaments are cosmic ray dominated (MWCCR and MWCCRd28); the mass fraction /m, however, 
is smaller. This is because of the smaller gas density and thermal pressure in the cold filaments. 
In addition, Figure [8] shows (short dashed line) that the filaments are more anisotropic {L\\/L± 
is larger) in the nonlinear phase for cosmic-ray dominated filaments; this is because the cosmic 
ray pressure resists parallel compression. Indeed, a visual comparison of the filaments with and 
without a large cosmic-ray pressure in Figure [12] shows that the absolute parallel length-scale of 
the filaments is larger when the cosmic rays are dynamically important. 

For a realistic cooling function, the density contrast between the filaments and the thermal 
plasma is much larger than in our simulations (because the stable thermal phase has T ~ 10"^ K 
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rather than T ~ 2 x 10^ K). For an initial ICM temperature of 10^ K, the real density contrast 
should be ~ 10'^ at a fixed pressure (assuming the filaments are not cosmic ray pressure dominated). 
Thus, even with an initially very small cosmic ray pressure in the ICM of pcr/p ~ 10~^, the cosmic 
rays can be adiabatically compressed to be dynamically important in filaments. This suggests 
that the cosmic ray dominated results in the left panel of Figure [12] are likely to be the most 
physically realistic. However, for the large gas densities and cosmic ray pressures that obtain in 
the filaments, cosmic ray losses due to ionization, pion production, and cosmic ray streaming will 
become important. The hadronic and ionization loss timescales ar e comparable, 2 00/nf.(cm~'^) 



Myr, for relativistic protons with kinetic energy of a few GeV (e.g.. ISchlickeiser II2002I ). The energy 
loss timescale because of cosmic ray streaming is roughly the Alfven crossing time along the filament 
(~ 1 Gyr for a 10 kpc long filament and an Alfven speed of 10 km s~^). Since these loss timescales 
are only modestly longer than the nominal cooling time, and since the filaments are expected to 
be dense, cosmic ray losses have to be included self-consistently. While including ionization and 
had ronic losses is straightforward, nu merically implementing cosmic ray streaming is non-trivial 



(see ISharma. Colella. Martin 1 12009| ) . A self-consistent treatment of this physics is beyond the 



scope of the present paper, but may modify the impact of cosmic rays on filament formation. 

The only non-adiabatic cosmic ray physics in our calculations is diffusion along magnetic field 
lines (Equation ([5])). Our calculations with different parallel diffusivities Du show that, so long as 



D\\ < 10^^ cm^s ^, the adiabatic results in Figure [T2] are reasonably applicable. ISharma et al. 



(j2009al ) presented general arguments that the diffusivity is likely to satisfy this inequality, so we 



suspect that large cosmic-ray pressures in filaments are the norm. There are indee d observational 



indic ations that this is the case (e.g., the modeling of atomic and molecular lines by lFerland et al 



20091 ): we will discuss this comparison in 



4.5. Simulations with Different Cooling/Heating Functions 

To understand which aspects of the nonlinear evolution of the thermal instability in the ICM 
are robust, we have carried out similar calculations to those reported here with different assumptions 
about the heating and/or cooling functions. Recall that the heating function is particularly poorly 
constrained in the ICM. To give one example, we carried out a series of simulations with the heating 
proportional to density, i.e., with a heating that is constant per unit mass instead of constant per 
unit volume as in our fiducial models shown here. The results were qualitatively similar to the 
fiducial case, with anisotropic filaments and most of the volume in the hot phase. Nonlinearly, the 
mass fraction in the cold phase (0.38 at 1.9 Gyr) is smaller than in the fiducial run. The cold phase 
is slightly hotter (?a 2 x 10^ K instead of ~ 10^ K in the case of the fiducial run; see Figure [6]) 
and has a smaller spread in temperature than in the fiducial run, and it takes longer for nonlinear 
saturation because the cold phase is heated more effectively than in the fiducial case. The aspect 
ratio of the cold filaments (measured by L\\/L±) is similar to the fiducial run. 

We also carried out simulations in which the volume averaged instantaneous heating rate was 
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not equal to the cooling rate: random perturbations (up to 200%) in both space and time were 
added to the volume averaged heating rate. These runs also showed results qualitatively similar 
to the fiducial run, except that the field lines were more disturbed from the initial configuration, 
and the mass and volume fractions of plasma at intermediate temperatures was larger (as would 
be expected). This demonstrates that the existence of a multiphase medium and cold filaments 
aligned along the magnetic field are robust consequences of thermal instability in the ICM. The 
only way out of these conclusions is if there is a heating mechanism that is locally thermally stable 
on scales > the Field length; this is a muc h more stringent r equirement to satisfy than the global 



thermal stability of the ICM (however, see lKunz et al. Il20ld ) 



The fiducial run (and all other runs) uses the modified cooling curve shown in Figure [T] with 
the stable phase at T < 2 x 10^ K. To assess what happens to filaments with a realistic cooling 
function, in which the stable phase is at < 10^ K, we carried out two runs (with and without 
cosmic rays) in which the stable phase of the cooling curve exists for T < 10^ K. The Field length 
in the stable phase of these simulations is ~ 8 times smaller than in our fiducial calculations (see 
Equation (fT3]) ). For this reason these runs are only barely resolved (see ^ 14.31 for discussion 
of convergence), but they nonetheless indicate the trends expected for a more realistic cooling 
function. Nonlinearly, the run with cosmic rays shows much longer (and broader) filaments than 
the run without cosmic rays. This can also be seen by comparing the filaments in Figure H] (for the 
fiducial run) and Figure [12] (for the run with cosmic rays MWCCR; see also Figure [8]); however, 
the difference is even more dramatic for the runs with a cooler stable phase. For a smaller stable 
phase temperature, simulations without cosmic rays have very narrow and short filaments, while 
in simulations with cosmic rays the sizes of filaments do not depend significantly on the stable 
phase temperature. Thus, adiabatically compressed cosmic rays, which dominate the pressure in 
the filaments, are likely able to prevent compression of the cold plasma to very small scales for a 
realistic cooling function . 



4.6. Runs v^ith Heating ^ Cooling 

In cluster cores the instantaneous heating rate is probably not identically equal to the cooling 
rate, as we have assumed in our models. However, the inferred global stability of clusters suggests 
that for timescales longer than a few cooling times heating does roughly balance cooling. Otherwise, 
all of the plasma will be in the cold phase (if cooling dominates) or in the hot phase (if heating 
dominates). To test the sensitivity of the phase structure to the degree of imbalance between heating 
and cooling, we carried out simulations in which heating does not quite balance cooling (Table [I]). 
Figure [13] shows temperature contour plots after ~ 0.95 Gyr for simulations with a constant heating 
per unit volume = 0.9xcooling (MWChO.Dc) and with heating per unit volume =1.05xcooling 
(MWChl.05c). The two plots differ dramatically. When cooling is somewhat stronger than heating 
(MWChO.Oc), the filaments are longer, much broader, and contain more of the mass, as compared 
to the fiducial run; by contrast, when heating exceeds cooling (MWChl.05c), the cold structures are 
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much smaller. The results differ even more dramatically from our fiducial calculations for a larger 
imbalance between heating and cooling. When heating does not exactly balance cooling, there 
will only be a single phase if we wait long enough. The relevant timescale is isec ^ icoo\C /\C — 
if I, the timescale for secular heating/cooling of the plasma, where H/C is the volume averaged 
heating/cooling rate. After ~ tgec the plasma will be dominated by hot/cold phase if heating/cooling 
dominates. Note that Figure [13] is shown at ~ 0.95 Gyr, which is ~ tscc for these models. The 
fact that many cluster cores show a multiphase structure implies that heating balances cooling 
over a few cooling times. In the future, a more quantitative comparison between simulations like 
those reported here and observations might provide interesting constraints on the degree of thermal 
balance in cluster cores. 



5. Astrophysical Implications 

Early thermal stability analyses of galaxy clusters were done within the context of the cooling 
flow model, in which mass inflows on the same times cale that the plasma co ols; this can signiflcantly 



modify the physics of the thermal instability (e.g.. iBalbus &: Soker Ill989l ). However, observations 
now clearly demonstrate the lack of signiflcant cooling flows; a poorly understood source of heating 
(plausibly a central AGN) roughly balances cooling, maintaining approximate global thermal sta- 
bility. In spite of their global stability, clusters can still be vulnerable to local thermal instability 
whenever Field's criterion is satisfled (i.e., whenever there are growing solutions to Equation dH])). 
Thermal conduction helps stabilize cluster plasmas on scales smaller than the Field length (Equa- 
tion ([9])), so it is the larger scale perturbations that are particularly prone to instability. It is not 
guaranteed that such instabilities in fact exist: whether they do depends on the details of how the 
plasma is heated. In this paper we have used two-dimensional MHD simulations with anisotropic 
conduction and cosmic rays to study the nonlinear dynamics of thermal instability for conditions 
appropriate to galaxy clusters, under the assumption that local heating is not able to maintain 
thermal stability. Our results can only be semi-quantitatively applied to observed clusters, given 
current uncertainties in the heating physics. Nonetheless, we find that none of our conclusions are 
that sensitive to the precise form of the heating function (e.g., whether it is constant per unit mass 
or constant per unit volume; ^4.5p . We also find similar results in simulations that include a slight 
imbalance between heating and cooling (so long as the simulation is not run too long; N4.6P or ran- 
dom perturbations in the heating/cooling rates on top of a thermal balance f ^4.5D. Observations 



of at o mic and molecular filaments and star formation in cool cluster cores (e.g., ICavagnolo et al. 



2008 



Q'Dea et al. II2008I ) provide observational evidence for local thermally unstable regions in 



clusters. 

Our calculations show that, for numerical convergence, the Field length in the cold medium 
needs to be resolved not only along the magnetic field, but also perpendicular to the field lines. To 
do so, we have artificially increased the temperature at which the plasma is thermally stable on the 
low temperature part of the cooling curve (to 2 x 10^ K; see Figured]). We have also included a 



- 20 - 



small isotropic thermal diffusivity, to ensure that the perpendicular Field length is resolved ( ^2.2p . 

During the evolution of the thermal instability, rapid thermal conduction along magnetic field 
lines suppresses compression of plasma along the field at scales smaller than the Field length. 
However, compression occurs perpendicular to field lines on large scales, where magnetic tension 
is not important. Thus if the Field length is < the size of a galaxy cluster core, and if the cooling 
time is short compared to the age of the cluster, the ICM is likely to be multiphase, with atomic 
filaments aligned with the local magnetic field. Note that this conclusion holds even in the fully 
nonlinear regime and does not require dynamically strong magnetic fields (Figure H]). Rather, 
thermal instability leads to a filamentary structure because of the poor heat transport across 
magnetic fields. This result implies that the orientation of atomic filaments can provide a local 
measure of the magnetic field direction in clusters. It also provides a physical ex planation for the 



filamentary structures se en in optical emission line observations of cluster cores (IConselice et al. 



2OOII : ISparks et al. II2004I ). Note that simulations with isotropic conduction show no preference for 



the cold gas to align with the magnetic field direction (e.g., Figure [7]). 

The filamentary structure in the cold gas is also imprinted on the diffuse X-ray emitting plasma 
in the hot ICM (e.g., FigureH]). Because of the large conductivity of the hot plasma (Equation (jlip ). 
it is natural for a given magnetic field line to become relatively isothermal. If different magnetic 
field lines undergo slightly different heating/cooling, as must surely be the case to some extent, this 
will lead to different temperatures, densities, and X-ray emissivities along different magnetic field 
lines. This could potentially expl ain the long, soft X-ray emitting isothermal structures observed 



in some clusters (ISun et al. II2009I ). 



The ambient cluster magnetic field is enhanced by flux freezing during the formation of fila- 
ments. Moreover, this enhancement of B extends over a region that is much longer than the extent 
of the cold gas itself (compare Figure [5] for \B\ with FigureH] for T). This is because as the plasma 
compresses along the magnetic field, it leaves behind regions devoid of much cool plasma that have 
nonetheless had field amplification by flux freezing. The volume averaged velocities induced by 
the thermal instability are ~ 25 km s~^ for our typical cluster parameters. The velocity in the 
hot phase is, on an average, directed toward the cold filaments; velocities in the cold filaments are 
larger (~ 100 km s~^), are generally parallel to the filaments, and may have strong shear (because 
of the merger of oppositely moving filaments; see the right panel of Figure [5]). The velocities ~ 100 
km s~^ we find in c old filaments are si milar to the measured random velocities of optical emission 



line filaments (e.g.. lHatch et al. 1120061 ). However, this comparison may be misleading because we 



ignore gravity in our simulations which can easily induce large motions in the dense filaments. 

The length of cold filaments along the magnetic field roughly scales with the Field length in 
the cold phase. The width of filaments is also determined by the perpendicular Field length. For 
a realistic atomic cooling curve in which the cold atomic gas is at ~ 10^ K, the filaments are 
expected to be extremely small ~ 10~^ pc. However, the observed atomic (e.g.. Ha) filaments are 
much longer than this. This can be explained if the filaments are supported by cosmic ray pressure 
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which prevents the collapse of the cold gas (see Figures [8] & [T2|) . The presence of a significant 
population of cosmic rays is also inferred by modeling the atomic and molecular lines from clusters 
(jFerland et al. II2009I ). Even if the cosmic ray pressure is small in the diffuse ICM {pcr/p ^ 10~^), 
adiabatic compression will result in cosmic ray dominated cold filaments. For the scales of interest, 
cosmic ray diffusion can be neglected if the c os mic ray diffusio i i coeffi cient is equal to the Galactic 



value (10^^ cm^s"^ iBerezinskn et al. Ill99d ): ISharma et al. I (|2009al ) argue that this is likely to 



be the case. Other loss processes (e.g., pion production, ionization, Alfven-wave excitation) may, 
however, be important, and could modify how effectively cosmic rays can support filaments ( ^4.4p : 
this will be studied in more detail in future work. If cosmic ray pressure is indeed substantial, the 
pressure of the thermal plasma in cold filaments can be significantly smaller than that of the ambient 
ICM. Substantial cosmic ray (and magnetic) pressure could i n principle help exp l ain the lack of 



star formation in the molecular filaments of NGC 1275 (e.g., iFabian et al. II2008I : I Salome et al. 



20061 ). Hadronic interactions between cosmic rays and thermal nucleons in dense filaments can 
produce a significant gamma-ray fiux due to neutral pion decay; however, because of the small 
volume occupied by the filaments, it is unlikely that the filaments will be detectable by current 
instruments. 

The perpendicular thickness of the filaments in our calculations is set largely by the isotropic 
diffusivity we include to ensure numerical convergence; in reality, however, the perpendicular ther- 
mal conductivity is negligible and some other physics (perhaps cosmic ray pressure again) must 
determine the perpendicular scale of the filaments. The properties of observed filaments can in 
principle be tested by Faraday rotation measured along and across the filaments; the Faraday rota- 
tion should be substantial (~ 10^ rad m~^ for rig ~ 10 cm~^, B ~ 10/iG, and a filament length of 
10 kpc). However, these observations are difficult precisely because the filaments are narrow and 
because this requires a reasonably strong radio source behind the filament. 

We have not included gravity in our simulations in order to focus on the thermal phase structure 
of the ICM and not the ICM dynamics. In stratified plasmas there will be a complex interplay 
between the thermal instability a r id buoyant motions (e.g., driven by buoyancy instabilities; see 
Parrish. Stone. &: Lemaster 1 120081 : iParrish. Quataert. &: Sharma I l2009l ) : this will be the focus of 
future work. In the nonlinear limit, cold, dense filaments are expected to fall, almost at the free 
fall rate, toward the cluster center. Magi ietic anchoring and levitation by unde r dense, buoyant 



bubbles may, however, prevent this (e.g., iHatch et al. 1120061 : iRevaz et al. M2008I ). Even with a 



significant gravitational field, we expect the filaments to be aligned with the local magnetic field 
as a consequence of the basic thermal physics of the ICM (i.e., cooling and anisotropic thermal 
conduction) . 
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A. Linear Stability Analysis 

We assume a background hydrostatic and thermal equilibrium. Let the equilibrium quantities 
{po, -Bo, Pcr,o) be constant in space; the following analysis is valid for kH » 1 where k is the 
wavenumber and H is the scale over which equilibrium quantities vary. We do not include gravity 
in the following analysis, and hence because of kH ^ 1 all terms involving background gradients 
are small. 

Perturbations of the form gi-™t+'tk-x) ^^^^ assumed, where w is the frequency. Linear pertur- 
bations are preceded by a 5 and equilibrium quantities have a subscript 0. The linearized equations 
are given by 

^ + ik-^ = 0, (Al) 
Po 

where ^ = iSv / uj is the displacement vector, 

- = _,fe^(^ + Pc. + i?V8vr) ^ ^^^^ ^^^^ 
Po 47rpo 
SB = i{k ■ Bo)$, - i{k ■ $)Bo, (A3) 

_ (tlZlM = -5[neniA{T)] - ik ■ SQ, (A4) 
7-1 

where = p/p, and the perturbation of space-constant H{t) vanishes as it equals volume averaged 
cooling rate which is constant in time in the linear regime, 

- iw{5pcr - IcrV^ crSp) = -ik ■ ST, (A5) 

where vh^ = Pcv/p- 



A.l. Fast Sonic Speed Limit 

Dotting Equation ()A2p with k gives 6pt/pt ~ (tsnd/^cooO^'^p/po; where Pt = P + Pcv + B'^ /Sir, 
*snd ~ ^{Pt/ pY^"^ ■, and w ~ i^ooi (^■^•' considering the condensation mode which grows at 
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the cooling time; t^^j = (7 — l)nenjA/po)- In the hmit when the sound-crossing time is shorter 
than the coohng time, the relative perturbation in total pressure is much smaller than the relative 
perturbation in density. Thus we can combine Equations ()A3p . ()A4p . and ()A5P to give, 

- iw6pt = -iw (7^1 + 7crf t,cr) " '^'^^ ( 8~ ) ~ ~ ^)^i'^eniA{T)] - (7 - l)ik SQ-ik ST. (A6) 



Now in the limit of isnd ^ ^cooi, the left hand side of Equation (|A6P can be ignored with respect to 
the first term on the right hand side. In this limit, from Equation (|A3p . we get 

diB^Sn) = vl6p/il - klv\/w^), (A7) 

where v\ = Bq/Attpq and /cy = fe • bg. Thus in the short sound-crossing time limit. Equation ()A6p 
reduces to 



5p = -(7 - l)S[neniA{T)] - (7 - l)ik ■ SQ - ik ■ ST, (A8) 



1 — /cyf^/tt;^ 

where pt is held constant in evaluating the right hand side. The cooling term can be written as 

(A9) 



d\neniA(T)], , d\neniA(T)] , „2 n 

5[nen,A{T)]p^ = ^^^-^\p6p - ^ ' ^ ^^ Ip^iPcr + B^Stt), 



dp dp 

where d[neniA{T)]/dp\p = {ueTiiT /p)dA/ dT , and d[neniA{T)]/ dp\p = -(neUiT'^ / p)d[A/T'^]/ dT . 
Thus, Equation (jA9|l . on combining with Equation (IA7I) . becomes 



,51n[neniA(r)] 



d In A 5pcr 
dlnT pq 



(iln(A/r2) 
dlnT 



+ 



dlnA 



P{l-klv\/w'^)d\nT 



5p_ 
Po' 



(AlO) 



where /3 = Sirp/B'^. From Equation ([7]) we obtain ifc • ST = D^^k'^^dpcr, so combining with Equation 
i5\i. we get 



Spc 



-iwjcrVlcr^p 

{-iw + ■ 



(All) 



Similarly, ik ■ SQ = Kk?,5T, and on using Equations ^A^ and (IATO]) . gives 



/ • ,2n<^?' -1 rflnAJpcr / . ,2-1 

(-zu;+X|| A^ii ) — = *cooi^]^ — + (-^^7 + X\\ ^11 + icooi 



(iln(A/r2) 



+ 



dlnA 



dlnT P{1 - kf.vyw^) dlnT 



6p 



(A12) 



where x\\ = {1 ~ is the thermal diffusivity. Combining Equations ()A12p . (jAlip . and (|A7 

and using 5pt ~ gives. 



-iwa 



+ t 



{-iw + D\\kp 
dln(A/r2 



cool 



-iw + x\\kl + tj, 
2 



dlnA^^ 2 i-iw + X\\kl) 
' ;5(1 - A:2t,2/«;2) 



™°i(ilnr 
dlnA 



-iwj + Xll^n 



(ihiT 



/3(1 - A:2t,2/?i;2) dliiT 



(A13) 
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where a = Pcr/p- In the hydro hmit (a ^ 1 and /3 ^ 1, irrespective of k'^\v\t~^^-^; i.e., magnetic 
tension plays no role in the condensation mode), we recover the classic isobaric thermal instability 
stabilized by conduction along field lines (Equation ([8|)). For /3 ^ 1 or a ^ 1 one obtains the 
thermal instability in the isochoric limit (see next section), with conductive stabilization for scales 
smaller than the isochoric Field length (Equation (jlOp ). The condensation mode is isochoric when 
magnetic/cosmic ray pressure dominates because the constancy of the total pressure is equivalent to 
the constancy of the magnetic/cosmic ray pressure, and from Equations (1A7P and (jAlip a constant 
magnetic/cosmic ray pressure implies a constant density. 

A. 2. Slaw Sonic Speed Limit 

In the opposite limit, tcooi ^ *snd) ^pI P ^ <^Pt/pt> Equation ()A6P gives 

- iwbpt = -(7 - l)5[neniA(r)] - (7 - \)ik ■ SQ - ik ■ ST, (A14) 

where terms on the right hand side are evaluated keeping the density constant (i.e., isochoric). The 
perturbed magnetic and cosmic ray pressure vanish in the isochoric limit (Equations ()A7p . (jAlip ) 
and 6pt = Sp, and the dispersion relation is the same as Equation (fTUj) . 
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Fig. 7. — Contour plots of LogiQ temperature (keV) for the simulation with isotropic thermal 
conduction at the Spitzer value (MWIC), at 0.475 Gyr (left) and 0.95 Gyr {right). The arrows 
show the magnetic field direction. 
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Fig. 8. — The anisotropy of the density field L||/L_|_ as a function of time for different runs: 
the fiducial run (MWC), the run with initial Pcr/p = 0.1 (MWCCR), and the run with isotropic 
conduction (MWIC). Note also that the filaments are longer and broader for simulations that 
include cosmic rays, i.e., both Ly & L±_ are larger even though L^jL^ is comparable (see Figures 

a&ii]). 




Fig. 9. — Contour plots of Logio temperature (keV) at 0.95 Gyr for higher (MWCh; left) and lower 
(MWCl; right) resolution analogues of our fiducial simulation. Figure H] shows the corresponding 
temperature plot for the fiducial run. All three are reasonably similar. 
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Fig. 10. — Contour plots of Logio temperature (keV) at 0.95 Gyr for high (2048; left) and low 
(512; right) resolution simulations without the small isotropic conductivity which is needed for 
convergence. Compare with Figure [9] which shows results for simulations including a small isotropic 
conductivity. 
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Fig. 11. — Left: Volume averaged filament anisotropy Equations & (fT5]) ) as a function 

of time. Right: Volume averaged random velocity as a function of time. In both cases, we show 
simulations with (labeled "Y") and without (labeled "N") a small isotropic conductivity (see i l2.2p . 
Simulations with the isotropic conduction converge reasonably well with increasing resolution (for 
A'" > 512; also see Figure [3|) but those without it do not. We could not run the higher (2048) 
resolution simulations for long because of very limiting time-step constraints. 
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Fig. 12. — Contour plots of the cosmic ray to plasma pressure ratio, LogiQ{pcr/p), at 0.95 Gyr for 
the runs with initial p^r/p = 0.1 (MWCCR; left) and initial pcr/p = 10"^ (MWCCRs; right). The 
density/temperature contour plots look similar to these because Pcr/p^^'^ is conserved. 
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Fig. 13. — Contour plots of LogiQ temperature (keV) at 0.95 Gyr for the run with volume averaged 
heating = 0.9 x volume averaged cooling (MWCh0.9c; left) and for the run with volume averaged 
heating = l.OSxvolume averaged cooling (MWChl.OSc; right). 



